Introdução à Análise Estatística Espacial

O que é Análise Estatística Espacial ?

  • São métodos estatísticos que levam em consideração a localização espacial do fenômeno estudado;

  • Segundo Bailey & Gatrell (1995), “Análise estatística espacial pode ser realizada quando os dados são espacialmente localizados e se considera explicitamente a possível importância de seu arranjo espacial na análise ou interpretação dos resultados”

  • A primeira pergunta a ser feita é: A distribuição dos dados apresenta um padrão aleatório ou apresenta uma agregação definida (clusters) ?

Origem da Estatística Espacial

  • Dr. John Snow (1813-1858) \(\rightarrow\) Considerado pai da Epidemiologia Moderna

Mapeamento dos casos de coléra (\(\bullet\)) e as bombas de água (X) em Londres, 1854.

Objetivos da Estatística Espacial

  1. Análise de padrões espaciais e espaço-temporais (ex: Análise Exploratória Espacial, Correlação Espacial)

  2. Modelagem Espacial (ex: Modelos Estatísticos de Regressão Espacial)

Dependência Espacial ou Autocorrelação Espacial

  • “Independência é um pressuposto muito conveniente que faz grande parte da teoria estatística matemática tratável. Entretanto, modelos que envolvem dependência estatística são frequentemente mais realísticos. . . . Neste tipos de dados a dependência está presente em todas as direções e fica mais fraca a medida em que aumenta a dispersão na localização dos dados.” (Cressie N.,1991)

  • “Todas as coisas se parecem, coisas mais próximas são mais parecidas que aquelas mais distantes.” (Tobler, 1979). Também conhecida como \(1^a\) Lei da Geografia

Tipologia dos Dados Espaciais

Os diferentes tipos de dados espaciais são tradicionalmente classificados de acordo com uma tipologia. Esta caracterização diz respeito a natureza estocástica da observação.

Noel Cressie (1993) divide a estatı́stica espacial em 3 grandes áreas:

  1. Dados de processos pontuais

  2. Dados de geoestatística

  3. Dados de área

Existem métodos estatı́sticos diferentes para descrever ou analisar estes tipos de dados. Exemplo:

Padrões Pontuais

  • O principal interesse está no conjunto de coordenadas geográficas representando as localizações exatas de eventos.

  • Exemplos: Localização de crimes, localização da residência dos casos de dengue, localização de espécies vegetais, etc.

  • Neste caso, o dado aleatório de interesse é a localização espacial do evento.

Estimativas de Kernel (ou Mapas de Calor)

  • Estimador de intensidade de distribuição de pontos:

\[\hat{\lambda}_{\tau}(u) = \dfrac{1}{\tau^2}\sum k(\dfrac{d(u_i , u)}{\tau}), d(u_i , u) \leq \tau\] *Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004

  • Localização da ocorrência de casos de dengue em Belo Horizonte/MG

Geoestatística

  • Podemos definir como sendo uma análise de um atributo espacialmente contínuo amostrado em localizações fixas.

  • Os dados compreendem um conjunto de localizações (em geral latitudes e longitudes), mas atribuidos a eles uma medida, como por exemplo: volume de chuva, concentração de poluentes no ar, número de ovos de Aedes postado em ovitrampas, etc.

  • A determinação experimental do semivariograma, para cada valor de \(x_i\), considera todos os pares de amostras \(z(x)\) e \(z(x+h)\), separadas pelo vetor distância \(h\), a partir da equação:

\[\hat{\gamma}(h) = \dfrac{1}{2N(h)}\sum^{N(h)}_{i=1}[z(x_i) - z(x_i + h)]^2 \]

  • Sendo \(\hat{\gamma}(h)\) o semivariograma estimado e \(N(h)\) o número de pares de valores medidos, \(z(x)\) e \(z(x+h)\), separados pelo vetor \(h\).

  • Esta fórmula, entretanto, não é robusta. Podem existir situações em que variabilidade local não é constante e se modifica ao longo da área de estudo (heteroscedasticidade).

*Fonte: Referência científica: Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds) “Análise Espacial de Dados Geográficos”. Brasília, EMBRAPA, 2004

  • Exemplo: Mapa sobre o teor de argila no solo.

Dados de Área

  • Na análise de áreas o atributo estudado é em geral resultado de uma contagem ou uma medida de síntese (em geral somas ou médias) em uma determinada área bem definida, por exemplo um polígono.

  • Tal polígono cuja forma pode ser complexa bem como as relações de vizinhança.

  • O objetivo é a detecção e explicação de padrões e tendências observados entre áreas.

Mapa Temático

  • Taxas de câncer de pulmão na população branca masculina nos Estados Unidos, por condados no ano de 1998

Matriz de Vizinhança

Moran Global, Moran Local e Lisa Map

Moran Global

Moran Global

Moran Local

Moran Local

  • Desigualdade no nível distrital na cobertura de saúde reprodutiva, materna, neonatal e infantil na Índia

Fonte: PANDA, Basant Kumar; KUMAR, Gulshan; AWASTHI, Ashish. District level inequality in reproductive, maternal, neonatal and child health coverage in India. BMC public health, v. 20, n. 1, p. 1-10, 2020.

Modelos de Regressão Espacial

  • Hipótese de independência das observações em geral é Falsa \(\rightarrow\) Dependência Espacial

  • Efeitos Espaciais \(\rightarrow\) Se existir forte tendência ou correlação espacial, os resultados serão influenciados, apresentando associação estatística onde não existe (e vice-versa);

  • Como verificar ? \(\rightarrow\) Medir a autocorrelação espacial dos resíduos da regressão (Índice de Moran dos resíduos).

  • Autocorrelação espacial constatada ! E agora ? \(\rightarrow\) Utilizar Modelos de regressão que incorporam efeitos espaciais:

    1. Modelos Globais: utilizam um único parâmetro para capturar a estrutura de correlação espacial \(\rightarrow\) ex: Spatial Error Models (CAR)

\[y_i = \beta_0 + \sum^{p}_{k} \beta_k x_{ik} + \varepsilon_i \text{ , sendo } \varepsilon_i = \lambda W + \xi\]

  • Os efeitos da autocorrelação espacial são associados ao termo de erro. Sendo \(W\) a matriz de vizinhaça, \(\varepsilon\) a componente do erro com efeitos espaciais, \(λ\) é o coeficiente autoregressivo e \(\xi\) é a componente do erro com variância constante e não correlacionada.

  • A hipótese nula para a não existência de autocorrelação é que \(\lambda = 0\), ou seja, o termo de erro não é espacialmente correlacionado.

    1. Modelos Locais: parâmetros variam continuamente no espaço ex: Geographically Weighted Regression (GWR)

\[y_i = \beta_{0}(u_i, v_i) + \sum^{p}_{k} \beta_k (u_i, v_i) x_{ik} + \varepsilon_i \text{ , sendo } (u_i, v_i) \text{ as coordenadas geográficas}\] Fonte: ARDILLY, Pascal et al. Manuel d’analyse spatiale. 2018.

Serão feitas duas aplicações de Estatística Espacial utilização o pacote estatístico R:

Aplicação I: Construindo mapas a partir da biblioteca geobr

  • Bibliotecas que serão utilizadas:
library(ggplot2)
library(dplyr)
library(viridis)
library(geobr)
library(sf)
library(maptools)
library(leaflet)
library(ggspatial)
  • Para acessar os dados dos limites territoriais de todos os estados brasileiros é necessário utilizar a função read_state.
brasil.ufs <- read_state(code_state = "all", year = 2019, showProgress = FALSE)
  • Primeiro, vamos fazer um gráfico apenas com as geometrias.
ggplot(brasil.ufs) + geom_sf()

  • Para a construção de um mapa onde cada estado recebe uma cor de acordo com a sua região geogrática, procedemos da seguinte forma:
ggplot(brasil.ufs) + geom_sf(aes(fill = name_region)) + labs(fill = "Região")

  • Agora utilizaremos dados relativos a porcentagem de municípios com rede de esgoto de acordo com a unidade da federação (fonte dos dados ) para fazer um gráfico. Vamos associar esses dados a tabela de acordo com a variável State e padronizaremos a porcentagem para variar de 0 a 100.
acesso_san <- data.frame(code_state = c(12, 27, 16, 13, 29, 23, 53, 32, 52, 21, 51, 
    50, 31, 15, 25, 41, 26, 22, 33, 24, 43, 11, 14, 42, 35, 28, 17), com_rede = c(0.273, 
    0.412, 0.313, 0.177, 0.513, 0.696, 1, 0.974, 0.28, 0.065, 0.191, 0.449, 0.916, 
    0.063, 0.731, 0.421, 0.881, 0.045, 0.924, 0.353, 0.405, 0.096, 0.4, 0.352, 0.998, 
    0.347, 0.129))

ggplot(mutate(left_join(brasil.ufs, acesso_san, by = "code_state"), com_rede = 100 * 
    com_rede), aes(fill = com_rede), color = "black") + geom_sf() + scale_fill_viridis(name = "Municípios com rede de esgoto (%)", 
    direction = -1) + xlab("Longitude") + ylab("Latitude") + annotation_scale(location = "bl") + 
    annotation_north_arrow(location = "br") + theme_minimal()

  • Uma forma alteranativa de apresentar esses mesmos dados se dá pela apresentação de círculos com raios proporcionais a porcentgem de municípios com rede de esgoto no centroide de cada geometria (nesse caso, UF).
coord_pontos <- st_centroid(mutate(left_join(brasil.ufs, acesso_san, by = "code_state"), 
    com_rede = 100 * com_rede))

ggplot(brasil.ufs) + geom_sf() + geom_sf(data = coord_pontos, aes(size = com_rede), 
    col = "blue", alpha = 0.65, show.legend = "point") + scale_size_continuous(name = "Municípios com rede de esgoto (%)") + 
    xlab("Longitude") + ylab("Latitude") + annotation_scale(location = "bl") + annotation_north_arrow(location = "br") + 
    theme_minimal()

  • Uma alternativa interativa para trabalhar com mapas é com a utilização do pacote leaflet
addCircleMarkers(addTiles(leaflet(data.frame(st_coordinates(coord_pontos), com_rede = coord_pontos$com_rede, 
    UF = coord_pontos$name_state))), ~X, ~Y, label = ~as.character(paste0(UF, ": ", 
    com_rede, "%")), labelOptions = labelOptions(textsize = "13px"), radius = ~com_rede/10, 
    fillOpacity = 0.5)

Créditos ao Thiago Mendonça

Aplicação II: Análise exploratória utilizando os dados de dengue em Dourados/MS

OBS: Os dados e as malhas geográficas utilizadas nessa apresentação, estão disponíveis no seguinte endereço: (link)

  • Biliotecas do R que serão utilizadas
library(readr)
library(tidyverse)
library(sf)
library(maptools)
library(spatstat)
library(tmap)
  • Lendo a tabela da população por setor censitário e os shapes files do contorno e por setor censitário de Dourados/MS
unzip("dados/pop2010.zip", exdir = "dados")
pop2010 <- read_csv("dados/pop2010.csv")

unzip("malhas/Setor_UTM_SIRGAS.zip", exdir = "malhas")
setor.sf <- read_sf("malhas/Setor_UTM_SIRGAS.shp", crs = 31981)

unzip("malhas/contorno.zip", exdir = "malhas")
contorno.sf <- read_sf("malhas/contorno.shp", crs = 31981)

OBS: Um aspecto muito importante dos dados espaciais é o sistema de referência de coordenadas (CRS) que é usado. Por exemplo, uma localização de (140, 12) não é significativa se você sabe onde está a origem e se a coordenada x está a 140 metros, quilômetros ou talvez graus de distância dela (na direção x). (link)

  • Fazendo um join com as tabelas com os setores censitários + população
setor.sf <- setor.sf %>% mutate(idsetor = as.numeric(CD_GEOCODI)) %>% inner_join(pop2010, 
    by = "idsetor")
  • Lendo e plotando os casos de dengue georreferenciados em Dourados/MS
unzip("dados/dengue_dourados.zip", exdir = "dados")
casos <- read_csv("dados/dengue_dourados.csv")
casos.sf <- st_as_sf(casos, coords = c("X", "Y"), crs = 31981)
ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = casos.sf, 
    color = "red", size = 1) + theme_void()

  • Lendo e plotando os os pontos de coleta de lixo georreferenciados em Dourados/MS
unzip("dados/lixo_dourados.zip", exdir = "dados")
lixo <- read_csv2("dados/lixo_dourados.csv")
lixo.sf <- st_as_sf(lixo, coords = c("X", "Y"), crs = 31981)
ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo.sf, 
    color = "blue", size = 1) + theme_void()

Como podemos observar, existem alguns pontos de coleta de lixo fora do contorno de Dourados/MS

  • Uma forma de ficarmos só com os pontos dentro do polígono é utilizando o comando st_intersection.
lixo2.sf <- st_intersection(contorno.sf, lixo.sf)
ggplot(setor.sf) + geom_sf(fill = "white", color = "black") + geom_sf(data = lixo2.sf, 
    color = "blue", size = 1) + theme_void()

  • Utilizando as informações dos casos (pontos) + do lixo (ponto) + população de cada setor censitário (mapa temático)
ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = casos.sf, color = "red", 
    size = 0.7, aes(colour = "Caso"), show.legend = "point") + geom_sf(data = lixo2.sf, 
    color = "salmon", size = 1, aes(colour = "Lixo"), show.legend = "point") + scale_fill_distiller(palette = "PuBu", 
    direction = 1) + scale_colour_manual(values = c(Caso = "red", Lixo = "slamon")) + 
    theme_minimal()

  • Agora serão construidos buffers de 500m de distância ao redor de cada ponto de coleta de lixo. Isso é interessante para verificar se os casos de dengue ocorrem em um raio de até 500m de distância dos pontos de coleta de lixo. Ou seja, a pergunta é, será que a distância dos pontos de coeta de lixo influenciam a ocorrência do caso de dengue ?

Buffers: São polígonos que contornam um objeto a uma determinada distância. Sua principal função é materializar os conceitos de “perto” e “longe”.

lixo_buffer <- st_buffer(lixo2.sf, 500)

ggplot(setor.sf) + geom_sf(aes(fill = pop)) + geom_sf(data = lixo_buffer, color = "gray", 
    fill = "transparent", size = 0.4) + geom_sf(data = casos.sf, color = "red", size = 0.7) + 
    scale_fill_distiller(palette = "PuBu", direction = 1) + scale_colour_manual(values = c(Caso = "red", 
    Lixo = "slamon")) + theme_minimal()

  • Represntando os casos e o lixo de forma interativa
tm_shape(setor.sf) + tm_borders("black") + tm_shape(casos.sf) + tm_dots("red") + 
    tm_shape(lixo.sf) + tm_dots("green") + tm_shape(lixo_buffer) + tm_borders("blue") + 
    tmap_mode("view")
  • Convertendo o dado de pontos (padrão pontual) para dados de área
## conta casos por setor 
casos.sf$contador <- 1 
casos <- setor.sf |>  
  st_join(casos.sf) |> 
  filter(CLASSI_FIN == 1) %>%  ## seleciona somente os casos confirmados
  group_by(ID) |> 
  summarise(casos = sum(contador))

st_geometry(casos) <- NULL  ## remove atributos de geometria

## numero de depositos de Lixo por setor 
lixo.sf$contador <- 1 

lixo <- setor.sf |>  
  st_join(lixo.sf) |> 
  group_by(ID) |> 
  summarise(lixo = sum(contador))

st_geometry(lixo) <- NULL ## remove atributos de geometria

# Inserindo as contagens dos casos e de pontos de coleta de lixo no atributo com a geometria.
setor.sf <- setor.sf |> 
  left_join(lixo,by = 'ID') |> 
  left_join(casos,by = 'ID')
  • Plotando o mapa temático dos casos por setor censitário
plot(setor.sf["casos"])

  • Plotando o mapa temático dos pontos de coleta de lixo por setor censitário
plot(setor.sf["lixo"])

  • Calculando a taxa de incidência e plotando o mapa temático dos pontos de coleta de lixo por setor censitário
setor.sf$tx <- (setor.sf$casos/setor.sf$pop) * 1000
setor.sf$tx[is.na(setor.sf$tx)] <- 0  # Transformando os missings em zero

summary(setor.sf$tx)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   0.000   1.736   4.065   4.311  56.399
  • Plotando a distribuição da incidência em Dourados/MS
library(wesanderson)
pal <- wes_palette("Moonrise3", 20, type = "continuous")

ggplot(setor.sf) + geom_sf(aes(fill = tx), color = "black") + scale_fill_gradientn(colours = pal) + 
    ggtitle("Taxa de incidência de Dengue") + theme_void()

  • Kernel por atributos

Vamos plotar o kernel por atributos referente a taxa de incidência de dengue em Dourados/MS.

Primeiramente é necessário dissolver os polígonos em formato sf para obter o contorno. Nesse caso queremos preservar o atributo AREA

dourados.contorno <- st_union(setor.sf)
plot(dourados.contorno)

dourados.w <- as.owin(st_geometry(dourados.contorno))
  • Extraindo os centróides dos polígonos em Dourados/MS
centroides <- st_centroid(st_geometry(setor.sf))

# Transformando em os centróides em formato sp
centroides.sp <- as.data.frame(as_Spatial(centroides))
names(centroides.sp) <- c("X", "Y")

plot(centroides.sp)

  • Colocando os pontos no formato sp
centroides.ppp <- ppp(centroides.sp$X, centroides.sp$Y, dourados.w)

plot(centroides.ppp, pch = 19, cex = 0.5)

  • Fazendo o kernel por atributo da taxa de detecção
kernel.tx <- density(centroides.ppp, 500, weights = setor.sf$tx, scalekernel = TRUE)
plot(kernel.tx)

  • Construindo a matriz de vizinhança para verificar a autocorrelação espacial
library(spdep)
viz <- poly2nb(setor.sf)
viz
Neighbour list object:
Number of regions: 284 
Number of nonzero links: 1726 
Percentage nonzero weights: 2.139952 
Average number of links: 6.077465
  • Iremos precisar da coordenadas dos centróides
setor.sp <- as(setor.sf, "Spatial")  # convertendo em formato sp
coord <- coordinates(setor.sp)  # coordenadas dos centroidas dos poligonos de dourados
class(setor.sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
  • Verificando a malha de conectividade da vizinhança de Dourados/MS
viz.sf <- as(nb2lines(viz, coords = coord), "sf")
viz.sf <- st_set_crs(viz.sf, st_crs(setor.sf))

# Plota o grafo de conectividade por contiguidade
mapa.viz <- ggplot(setor.sf) + geom_sf(fill = "salmon", color = "white") + geom_sf(data = viz.sf) + 
    theme_minimal() + ggtitle("Vizinhança por \n conectividade") + ylab("Latitude") + 
    xlab("Longitude")
mapa.viz

  • Obtendo a correlação da taxa de incidência de dengue Dourados/MS
pesos.viz <- nb2listw(viz)
moran.test(setor.sf$tx, pesos.viz)

    Moran I test under randomisation

data:  setor.sf$tx  
weights: pesos.viz    

Moran I statistic standard deviate = 15.724, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.524720303      -0.003533569       0.001128660
  • Plotando o correlograma
correl <- sp.correlogram(viz, setor.sf$tx, order = 8, method = "I")
correl
Spatial correlogram for setor.sf$tx 
method: Moran's I
           estimate expectation    variance standard deviate Pr(I) two sided
1 (284)  0.52472030 -0.00353357  0.00112866          15.7239       < 2.2e-16
2 (284)  0.23776816 -0.00353357  0.00048540          10.9525       < 2.2e-16
3 (284)  0.09536593 -0.00353357  0.00028358           5.8729       4.282e-09
4 (284) -0.02063378 -0.00353357  0.00019736          -1.2172         0.22351
5 (284) -0.07589158 -0.00353357  0.00016732          -5.5939       2.221e-08
6 (284) -0.06448448 -0.00353357  0.00016165          -4.7940       1.635e-06
7 (284) -0.02708340 -0.00353357  0.00016725          -1.8210         0.06861
8 (284) -0.02744543 -0.00353357  0.00018328          -1.7662         0.07735
           
1 (284) ***
2 (284) ***
3 (284) ***
4 (284)    
5 (284) ***
6 (284) ***
7 (284) .  
8 (284) .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(correl)

  • Mapeando os polígonos que tiveram os p-valores mais significativos no Moran Local.
setor.sf$pval <- localmoran(setor.sf$tx, pesos.viz)[, 5]

tm_shape(setor.sf) + tm_polygons(col = "pval", title = "p-valores", breaks = c(0, 
    0.01, 0.05, 0.1, 1), border.col = "white", palette = "-Oranges") + tm_scale_bar(width = 0.15)
  • Moran Local (Lisa Map) da taxa de incidência Dourados/MS
resI <- localmoran.sad(lm(setor.sf$tx ~ 1), 1:length(viz), viz, style = "W")
summary(resI)[1:10, ]
    Local Morans I Stand. dev. (N)   Pr. (N) Saddlepoint  Pr. (Sad)
1 1   -0.009885292     -0.01575255 0.5062841 -0.03600714 0.51436167
2 2    0.226981101      0.46511580 0.3209243  0.70056077 0.24178858
3 3   -0.010910944     -0.01829621 0.5072987 -0.04041673 0.51611955
4 4    0.484675519      1.31014141 0.0950740  1.43930439 0.07503215
5 5    0.018648578      0.05952726 0.4762661  0.09384037 0.46261798
     Expectation  Variance    Skewness Kurtosis   Minimum  Maximum
1 1 -0.003533569 0.1625854 -0.05147857 8.753481 -57.75455 56.75455
2 2 -0.003533569 0.2456263 -0.04188294 8.752890 -70.87400 69.87400
3 3 -0.003533569 0.1625854 -0.05147857 8.753481 -57.75455 56.75455
4 4 -0.003533569 0.1388594 -0.05570264 8.753780 -53.41199 52.41199
5 5 -0.003533569 0.1388594 -0.05570264 8.753780 -53.41199 52.41199
            omega       sad.r       sad.u
1 1 -0.0001369880 -0.01569409 -0.01569909
2 2  0.0028119325  0.44472643  0.49831660
3 3 -0.0001590844 -0.01822753 -0.01823490
4 4  0.0066304485  1.08297547  1.59298211
5 5  0.0005595065  0.05929966  0.05942125
 [ reached 'max' / getOption("max.print") -- omitted 5 rows ]
setor.sf$MoranLocal <- summary(resI)[, 1]

library(scales)

ggplot(setor.sf) + geom_sf(aes(fill = MoranLocal), color = "black") + scale_fill_gradientn(colours = c("blue", 
    "white", "red"), values = rescale(c(min(setor.sf$MoranLocal), 0, max(setor.sf$MoranLocal))), 
    guide = "colorbar") + ggtitle("Moran local") + theme_void()

Bibliografia sugerida

  • Druck, S.; Carvalho, M.S.; Câmara, G.; Monteiro, A.V.M. (eds). Análise Espacial de Dados Geográficos. Brasília, EMBRAPA, 2004.

  • Interactive Spatial Data Analysis by Trevor C. Bailey , Anthony C. Gatrell Routledge, 1995

  • Applied Spatial Statistics for Public Health Data; Lance A. Waller, Carol A. Gotway Wiley-Interscience 1St ed. 2004

  • Applied Spatial Data Analysis with R; Roger S. Bivand, Edzer Pebesma , Virgilio Gomez-Rubio Springer; Edição: 2nd ed. 2013

  • Geocomputation with R by Robin Lovelace, Jakub Nowosad and Jannes Muenchow. Geocomputation with R, 2021.

Programa de Pós-Graduação em Epidemiologia em Saúde Pública (ENSP/FIOCRUZ)

  • Programa de Pós-Graduação em Epidemiologia em Saúde Pública (ENSP/FIOCRUZ) link)

Programa de Pós-Graduação em Ciências Veterinárias (IV/UFRRJ)

  • Programa de Pós-Graduação em Ciências Veterinárias (IV/UFRRJ) link

Obrigado !!!!!